HIV incidence and impact of interventions among female sex workers and their clients in the Middle East and north Africa: a modelling study

Summary Background The incidence of HIV infection among female sex workers and their clients in the Middle East and north Africa is not well known. We aimed to assess HIV incidence, the contribution of heterosexual sex work networks to these numbers, and the effect of interventions by use of mathematical modelling. Methods In this modelling study, we developed a novel, individual-based model to simulate HIV epidemic dynamics in heterosexual sex work networks. We applied this model to 12 countries in the Middle East and north Africa that had sufficient data to estimate incidence in 2020 and the impact of interventions by 2030 (Algeria, Bahrain, Djibouti, Iran, Libya, Morocco, Pakistan, Somalia, South Sudan, Sudan, Tunisia, and Yemen). Model-input parameters were provided through a systematic review of HIV prevalence, sexual and injecting behaviours, and risk group size estimates of female sex workers and clients. Model output was number of incident HIV infections under different modelling scenarios for each country. Summary statistics were generated on these model output scenarios. Findings Based on the output of our model, we estimated a total of 14 604 (95% uncertainty interval [UI] CI 7929–31 819) new HIV infections in the 12 countries in 2020 among female sex workers, clients, and spouses, which constituted 28·1% of 51 995 total new cases in all adults in these 12 countries combined. Model-estimated number of new infections in 2020 in the 12 countries combined was 3471 (95% UI 1295–10 308) in female sex workers, 6416 (3144–14 223) in clients, and 4717 (3490–7288) in client spouses. Contribution of incidence in heterosexual sex work networks to total incidence varied widely, ranging from 3·3% in Pakistan to 71·8% in South Sudan and 72·7% in Djibouti. Incidence in heterosexual sex work networks was distributed roughly equally among female sex workers, clients, and client spouses. Estimated incidence rates among female sex workers per 1000 person-years ranged from 0·4 (95% UI 0·0–7·1) in Yemen to 34·3 (17·2–59·6) in South Sudan. In countries where HIV acquisition through injecting drug use creates substantial exposure for female sex workers who inject drugs, estimated incidence rates per 1000 person-years ranged from 5·1 (95% UI 0·0–35·1) in Iran to 45·8 (0·0–428·6) in Pakistan. The model output predicted that any of the programmed interventions would substantially reduce incidence. Even when a subpopulation did not benefit directly from an intervention, it benefited indirectly through reduction in onward transmission, and indirect impact was often half as large as the direct impact. Interpretation Substantial HIV incidence occurs in heterosexual sex work networks across the Middle East and north Africa with client spouses being heavily affected, in addition to female sex workers and clients. Rapid scaling-up of comprehensive treatment and prevention services for female sex workers is urgently needed. Funding Qatar National Research Fund (a member of Qatar Foundation), the Biostatistics, Epidemiology, and Biomathematics Research Core at the Weill Cornell Medicine-Qatar, Qatar University-Marubeni, the UK Medical Research Council, and the UK Department for International Development.

: Map of the countries included in this study. Map includes the model-estimated HIV incidence * arising in heterosexual sex work networks for the year 2020 and its contribution to total HIV incidence in the population for 12 countries with sufficient data to conduct the simulations. .. 23 Table S2 Section S1. Details of the mathematical modelling methods

Heterosexual sex work network
In the model, each female sex worker (FSW) or client in the network enters/exits the sexual network, forms/dissolves sexual partnerships, or acquires HIV through sex or by injecting drugs at event-specific probabilities at each time step in each simulation run. The sexual network is constructed assuming that the number of sexual partnerships formed by each regular or nonregular client with FSWs follows a gamma distribution, reflecting sexual network and behavior studies. 1- 6 The mean and variance (at 25% of the mean) of these distributions were informed by country-level data on sexual behavior in heterosexual sex work networks (HSWNs) for each of regular and non-regular clients, with higher mean and variance for non-regular clients. 1 Each month, every regular or non-regular client may form a new partnership with one or more FSWs, based on a random probability drawn from these distributions. Existing partnerships may also dissolve stochastically assuming an exponential distribution at a rate of inverse of duration of partnerships, which varies based on whether they involve a regular or non-regular client.
Accordingly, in such sexual networks, each client randomly selects FSW partners, but clients may have different propensities to form partnerships, a situation known as proportionate mixing. 5,7 FSWs exit the HSWN if they cease to practice sex work, and for clients if they cease seeking sex with FSWs, or through natural and AIDS-related mortality (table S1). Lower HIV transmission, slower AIDS disease progression, and higher life expectancy were assumed for individuals on antiretroviral therapy (ART; table S1). Those who exit the HSWN are replaced by susceptible persons, thus maintaining a fixed cohort size for FSWs and clients. 6 While the model assumes that HIV acquisition among FSWs can occur through sex with a client or through injecting drug use with an injecting partner, HIV acquisition among clients was assumed to occur only through sex with an HIV-positive FSW. Other sources of infection, such as the client's spouse, other heterosexual partners, male same-sex partners, and injecting drug use were not considered. Evidence suggests that the risk of HIV infection through these modes of exposure among clients is probably substantially smaller than the risk of infection through sex with a FSW in most countries of the Middle East and North Africa (MENA). 1, [8][9][10]

HIV sexual transmission in FSW-client partnerships
Probability of HIV sexual transmission in an HIV sero-discordant partnership, that includes an HIV-positive FSW/client and a susceptible counterpart, was determined from the probability of transmission per coital act per HIV stage of infection, number of coital acts per partnership, which varied based on whether partnerships were with regular or non-regular clients, and interventions that affect HIV transmission.
These interventions included ART among FSWs and their injecting partners, condom use in the partnership, male circumcision in the client, and pre-exposure prophylaxis (PrEP) in the FSW.
Coverage of these interventions for FSWs and clients was based on data for each country and was implemented in the model by random assignment.

HIV transmission through drug injection
Proportions of FSWs who inject drugs were based on data for each country. Since the injecting partners of FSWs can be other FSWs or persons external to the modeled sexual network of FSWs and clients, HIV acquisition through injecting drug use was modeled through an external hazard rate (force of infection) that depended on whether the FSW was on PrEP and whether her 7 injecting partner was on ART. Otherwise, a constant hazard rate was assumed and was derived by fitting model output to country-level data on HIV prevalence among FSWs who inject drugs, 1 or alternatively if such data were not available, to HIV prevalence among people who inject drugs (PWID). 11 FSWs who inject were assumed to inject for a specific duration, set at 10 years, 11 which differed from the duration of sex work set at 35 years. 1 Effect of PrEP on injecting drug use transmission was assumed to be equivalent to that for sexual transmission.
The proportion of FSWs who inject drugs (table 2 of main text) is ≥10-fold higher than that among clients and spouses, assuming that the latter proportion can be approximated by that among the wider population. 11,12 Therefore, HIV acquisition through injecting drug use was not modeled among clients and spouses at it likely plays a much less significant role in transmission dynamics among them.

HIV sexual transmission from clients to their spouses
HIV sexual transmission from clients to their spouses was modeled using a separate deterministic model, but using the individual-based model output as input (section S2). This was done for computational efficiency and because the stochastic and non-linear effects are less likely to be critical here, as the transmission is one-directional from clients to their spouses.
Numbers of HIV transmissions from clients to spouses were estimated using the proportion of clients in spousal partnerships, HIV prevalence among clients, numbers of susceptible spouses, and probability of HIV transmission per partnership. The latter was estimated using the

HIV natural history
HIV natural history was based on established empirical epidemiological measures (table S1).
Progression through each of HIV infection stages was modeled assuming an exponential distribution through rates derived as the inverse of duration of each HIV stage and implemented through a stochastic process.

Data sources and model parameters
The primary data source for this modeling study was the recently completed systematic review of HIV, sexual and injecting behavior, and population size estimates in FSWs and clients in MENA. 1 Countries were included in the present study if they had sufficient input data to simulate the HIV epidemic in the HSWN and HIV prevalence among FSWs was ≥0.5%. Country-specific parameter values were selected based on the most recent representative studies identified through the aforementioned systematic review. 1 Priority was given to studies with rigorous sampling methodologies, such as integrated bio-behavioral surveillance (IBBS) surveys.
Where several nationally representative estimates based on IBBS surveys were available, 1 the 9 mean of these estimates was considered. Otherwise, data collected after the year 2000 were pooled using random-effects meta-analysis. This methodology used Freeman-Tukey type arcsine square-root transformation to stabilize variances 13,14 before weighting measures using the inverse-variance method, 14,15 followed by pooling using DerSimonian-Laird random-effects models to account for sampling variation and true heterogeneity. 16,17 Data for coverage of interventions were primarily based on findings of the systematic review, 1 or alternatively, on UNAIDS compilations, 18 or imputed using the regional median for these parameters. 1 Demographic and Health Survey data on men in the general population were used to derive, for each country, the proportion of clients in spousal partnerships (defined as a marital/cohabiting partnership for ≥1 year) and the proportion of sexual acts protected by condom use in these partnerships. 19 The proportion of FSWs in spousal partnerships was based on the extracted data of the systematic review of sexual behavior of FSWs. 1 For countries with missing information, measures were imputed by pooling regional data using random-effects meta-analysis.
The population size of FSWs and clients in each country was based on country-level data. 1 Other model parameters, such as for HIV transmission and efficacy/effectiveness and coverage of interventions, were based on current evidence in the literature (table S1 and tables 1-2 of main text).

Model simulations
The model-generated sexual network was established with a "burn-in" of 50 years to ensure equilibrium of network structure prior to HIV introduction. Subsequently, HIV infection was seeded and the model was run for an additional "burn-in" of 300 years to ensure epidemic equilibrium in each country by 2020. These long durations of burn-in were chosen to ensure the stability of the epidemic simulations. Shorter durations could have been chosen, but we opted for 10 caution to ensure no temporal effects bias our estimates. Since epidemiological measures of interest, such as HIV incidence, were estimated over a short time horizon of one year, and in absence of quality country-level trend data for HIV prevalence in FSWs and clients in nearly all MENA countries, 1 analyses were implemented starting from this epidemic equilibrium.
Model predictions for each country were based on the mean and 95% uncertainty intervals (UIs) of distributions of outcome measures generated by 500 simulation runs. At this number of runs, the mean and distribution of outcome measures minimally varied with an increase in the number of simulations. The 95% UIs were generated directly from these distributions of simulation runs after excluding runs with HIV stochastic extinction. Extinctions were excluded because they resulted from the finite size of the simulated sexual network (stochastic fade-out), but unlikely to occur in reality considering that the actual size of the sexual network is much larger. For computational efficiency, after experimenting with different cohort sizes, simulations were performed using a cohort of 600 FSWs and 6,000 clients (one-third of which are regular and two-thirds are non-regular/one-time clients), as informed by MENA data. 1 Outcome measures were subsequently scaled-up to reflect the actual population sizes in each country. 1

Model fitting
Model fitting to HIV prevalence data among FSWs and HIV prevalence among FSWs who inject drugs was performed to estimate the overall sexual partnership formation rate and the baseline hazard rate of acquiring HIV through injecting drug use in each included country. Nonlinear least-square fitting using the Nelder-Mead simplex algorithm 20 was implemented iteratively to generate a set of 50 best model fits. A best model fit was defined as a relative error of <5% between model predictions and empirical data. The final best model fit was the most probable value for the sexual partnership rate and injecting hazard rate among the 50 best model fits- FSWs, clients, and client spouses over the duration of a year, by the total HIV incidence in the population (15-49 years) during that year, as estimated by UNAIDS. 18

Impact of interventions
The impact of expanding HIV interventions among FSWs on HIV incidence arising in HSWNs was assessed by estimating, using 500 simulation runs, the mean number of infections that would be averted over a 10-year duration after implementing the interventions, and the proportional decrease in incidence during this time (table 3 of main text).

Oversight
Ethical approval was not needed for this modeling study because the study uses published data for model input and provides aggregate measures for HIV incidence and impact of interventions.

Section S2. Estimation of HIV incidence in stable partners (spouses) of clients of FSWs
We modelled HIV sexual transmission from clients of FSWs to their stable partners (spouses) HIV incidence rate is thus given by: ( parallel, as well as instances pre-configured with GPU libraries and common applications for computational analyses. Our sensitivity analyses indicated that a three-fold increase in the size of the simulated populations is optimal for reducing stochastic fade-out (table S2). However, generating only baseline estimates using 100 simulation runs for one country entailed running the model for 5-6 days thus undermining the feasibility of extending this increase in population sizes to all analyses.

Section S4. Limitations and caveats
Estimates were generated through mathematical modeling, and therefore should be seen as an approximation of actual reality based on our current knowledge and the available input data that parametrized the model. Analyses were possible for only 12 of 23 MENA countries with sufficient HIV prevalence, behavioral, and risk group size estimate data to apply the model. To minimize bias in input data, we systematically pooled these data for more robust mean estimates that account for the totality of available evidence. Some model input data were global rather than MENA-specific such as the real-world effectiveness in achieving viral suppression among FSWs. 38 Some possible effects were not included, such as sero-sorting or differences in uptake of interventions for different groups of FSWs and clients, as supporting evidence was not available to include them in the model. Effect of other sexually transmitted infections on HIV transmission probability was not also included due to conflicting evidence, 56  prevalence among FSWs is ≤0.5%. These countries were thus excluded from analysis (n=6). This may also have slightly underestimated HIV incidence in included countries due to finite-network effects and higher likelihood of stochastic extinction (stochastic fade-out). Sensitivity analyses were conducted where the sample sizes of the simulated cohorts were increased by 2-fold, 3-fold, and 4-fold, to investigate the effect of stochastic fade-out on simulated model outcomes. The stochastic fade-out had a minimal effect on the point estimates of the model outcomes, but the 95% UIs were narrower with the larger simulated cohorts (table S2). 22 The stochastic fade-out further resulted in higher stochasticity in simulations assessing the impact of interventions up to 2030. The impact was thus assessed after 30 years "burn-in" to reduce stochasticity, and then scaled back to a 10-year duration, which may have overestimated the indirect impact of interventions on onward transmission of infection. The indirect impact of interventions on incidence is slower to materialize than the direct impact. The latter, such as for condom use, is immediate the moment a condom is used in a simulated sexual partnership.
Both point estimates and uncertainty ranges were directly obtained from model simulations, and therefore may not have factored other sources of uncertainty that may exist in real world sexual networks. Uncertainty intervals around estimates of the relative contribution of HIV incidence arising in HSWNs to total HIV incidence were often too broad. This is due to the finite-network effects and stochastic fade-out, more so when the number of infections was small leading to a higher effect of stochasticity on predicted outcomes. Figure S1: Map of the countries included in this study. Map includes the model-estimated HIV incidence * arising in heterosexual sex work networks for the year 2020 and its contribution to total HIV incidence in the population for 12 countries with sufficient data to conduct the simulations.
* Incidence here refers to the number of new HIV infections in 2020.   27 Figure S4: Contribution of HIV incidence occurring in heterosexual sex work networks to total HIV incidence in MENA.  628 (53.1) Abbreviations: ART: antiretroviral therapy; FSWs: female sex workers; e: effectiveness; NA: not applicable; PrEP: pre-exposure prophylaxis; PWID: people who inject drugs. * Estimates for the number of averted infections have been rounded to the nearest digit and may not exactly match the corresponding proportion of averted infections. † Includes expanding ART coverage to 50% with efficacy in preventing HIV transmission to partners of 96%, increasing condom use to 50%, and increasing PrEP to 25%. Baseline coverage was used whenever it was higher than that set in the investigated scenario. ‡ Includes expanding interventions to the highest modelled coverage levels including expanding, ART coverage to 81% with efficacy of 96%, increasing condom use to 80%, and increasing PrEP to 50%.